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Several, recently proposed methods of surface manufacturing based on ion beam sputtering, which 
involve dual beam setups, sequential application of ion beams from different directions, or sample 
rotation, are studied with the method of kinetic Monte Carlo simulation of ion beam erosion and 
surface diffusion. In this work, we only consider erosion dominated situations. The results are 
discussed by comparing them to a number of theoretical propositions and to experimental findings. 
Two ion-beams aligned opposite to each other produce stationary, symmetric ripples. Two ion 
beams crossing at right angle will produce square patterns only, if they are exactly balanced. In 
all other cases of crossed beams, ripple patterns are created, and their orientations are shown to be 
predictable from linear continuum theory. In sequential ion beam sputtering we find a very rapid 
destruction of structures created from the previous beam direction after a rotation step, which leads 
to a transient decrease of overall roughness. Superpositions of patterns from several rotation steps 
are difficult to obtain, as they exist only in very short time windows. In setups with a single beam 
directed towards a rotating sample, we find a non-monotonic dependence of roughness on rotation 
frequency, with a very pronounced minimum appearing at the frequency scale set by the relaxation 
of prestructures observed in sequential ion beam setups. Furthermore we find that the logarithm of 
the height of structures decreases proportional to the inverse frequency. 

PACS numbers: 68.55.-a,68.35.Ct,05.45.-a,79.20.Rf,81.16.Rf 
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I. INTRODUCTION 

In recent years ion-beam sputtering (IBS) as a method 
of nanoscale surface manufacturing^ has attracted much 
experimental work (for a recent review see Ref.'^). Self- 
organized ripple patterns appear generically by oblique 
incidence single-beam IBS, and dot or hole patterns are 
obtainable by normal-incidence IBS or by oblique inci- 
dence IBS on rotating samples. To produce a larger vari- 
ety of structures and to improve their quality, more com- 
plex setups with multiple ion beams, IBS on prestruc- 
tured samples, and rotating samples have been used, but 
few of them have been investigated in detail. In partic- 
ular, Carter^^ has proposed the use of dual ion beam 
sputtering (DIBS), Vogel and Linz^ proposed a four- 
beam setup and claimed that corresponding results may 
be obtained from a sequence of prestructuring and step- 
wise beam or sample rotation using a single ion beam 
(sequential ion beam sputtering or SIBS). Continuous 
rotation of the sample or the ion beam (referred to as 
RIBS, i.e. rotating IBS, in the following) has been pro- 
posed to suppress ripple formation^^^ (for example in 
SIMS and AES, where ripple formation would reduce the 
depth profiling resolution), but also to enhance the qual- 
ity of isotropic structures.^ These proposals were based 
on continuum theories of the self-organized pattern for- 
mation by IBS, which are further developments of the 
work of Sigmund^^ , Bradley and Harper^ ^ and Makeev 
et al.jii^ Recently, Joe, Kim et al.^^ conducted a system- 
atic experimental study of DIBS and SIBS on Au(OOl). 
A number of their findings are not in accordance with 
expectations derived from the conventional continuum 
framework. This motivated us to study DIBS, SIBS and 



RIBS setups with Monte Carlo (MC) simulation meth- 
ods. Recent proposals^' are based upon extensions of 
the standard continuum model (which is formulated as an 
anisotropic and noisy Kuramoto-Sivashinsky equation). 
These extensions have introduced new physical mecha- 
nisms, which change the scenario of pattern formation of 
the standard model, but none of these mechanisms has 
been confirmed and tested independently beyond doubt. 
In this situation, it may be helpful to see, what can 
be achieved from a simulation, which is a straightfor- 
ward stochastic implementation of Sigmund's energy de- 
position and sputtering formula, combined with indepen- 
dently tested models of surface diffusion in simple model 
systems, but without any further approximations. 

In the next section, we briefly introduce our MC sim- 
ulation method and the geometries of the considered se- 
tups. Subsequently, we present and discuss results we 
obtained for the topographies, the shape and orientation 
of the ripples, the structure function and the evolution of 
roughness in different setups. We compare our findings 
to theoretical proposals, experimental results and to the 
standard linear continuum model. 



II. MONTE CARLO SIMULATION 

Throughout this work, we will focus on the erosion 
dominated regime, leaving more complex interplays be- 
tween erosion and diffusion in multi-beam and rotating 
setups for further studies. Parameters of the simulations 
(given below) are chosen appropriately, such that in sin- 
gle ion beam setups, ripples perpendicular to the direc- 
tion of the beam are produced on non-rotating samples 
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FIG. 1: Three different IBS experiment setups, (a) dual ion 
beam with fixed sample (DIBS), (b) continuously rotating 
sample with single ion beam (RIBS), (c) sequentially sputter- 
ing of pre-rippled surfaces after 90° rotation (SIBS). 

and cellular structures occur on rotating samples . -^^i-^^ 

Fig. [1] depicts the geometries of setups we use in the 
simulations. A DIBS setup consists of two ion beams 
incident from directions described by polar angles 6i and 
02 and azimuthal angles 0i and 02 (see Fig. [TJa)). As a 
simple special case, we will consider opposing beams, i.e. 
Oi = $2 and A(j) = (/>i — (/>2 = 180°. Crossed beams are 
studied for equal polar angles as well as for the general 
case of different polar and azimuthal angles. 

SIBS setups are shown in Fig. [D^c). In a first step 
ripples are produced by a single ion beam. Then we 
change the azimuthal angle of the ion beam direction by 
some Acj) and monitor the further evolution of surface 
structures. 

Finally, Fig.lTJb) shows the RIBS setup, which is char- 
acterized by a constant angular velocity uj of the sample, 
which we realize by a fixed sample and a correspondingly 
rotation beam. The evolution of structures and surface 
roughness will be systematically studied as a function of 

UJ. 

A first version of our Monte Carlo model of erosion and 
surface diffusion has been introduced in Ref. Il7|. Further 
developments of the model are discussed in Ref.[l8|. Here, 
we use a solid-on-solid cubic lattice with an initially flat 
surface of size L x L with L = 512. Its configuration 
is described by a single valued, time dependent, height 
function h{x^ t), which obeys periodic boundary condi- 
tions. 

Ions start from a random position above the surface 
and move towards the surface in direction (^, (j)). Ero- 
sion occurs as a result of energy transfer via collision 
cascades induced by impinging ions. The ions penetrate 
the surface through a depth of a. The average energy 
distribution of a cascade initialized by a single ion is ap- 
proximated by a Gaussian distribution which has width 
of a in directions of ion-beam and ja perpendicular to the 
ion-beam, respectively. As default values, we use a = 9, 



<j = 3 and /i = 1.5 lattice constants. Every lateral atom 
is a candidate for erosion, the probability of erosion be- 
ing proportional to the amount of energy, which reaches 
the atom. The parameters are tuned so that the average 
yield is ~ 7 atoms per ion. The natural time scale in this 
type of MC simulations is proportional to ion fluence, 
therefore we use 1 ML = ions, corresponding to 1 ion 
per atom at the flat surface as unit of time in this work. 
It is connected to laboratory time via the ion flux. 

Surface diffusion is modeled as thermally activated 
hopping. In a diffusion sweep every surface atom has the 
chance to hop to a nearest neighboring site. The attempt 
frequency is calibrated to a temperature of T AOOK 
and a substrate barrier of Es = 0.75eV is chosen. This 
implies that after every L^/1000 incoming ions, one dif- 
fusion sweep over the surface is taken. The hops are 
accepted with the probability of exp(— A^^/Zc^T). The 
additional energy barrier {AE) consists of two terms, 
= nEnn + EsB' Here n is the net number of bro- 
ken bonds after the hop. If this net number is negative, 
i.e. additional bonds have been created, n is put equal 
to zero. We choose Enn = O.lSeV as nearest neighbor 
bond energy. Esb = 0.15eV is an Ehrlich-Schwoebel 
barrier term (ES), which is only nonzero, if an ad-atom 
approaches a step edge from an upper terrace. More de- 
tails about the surface diffusion modeling can be found 
in Ref. 18. 

The rotation of sample during IBS is equivalent to a 
fixed sample and a rotating ion beam. Therefore rotation 
or any change of azimuthal angles in lab coordinates is 
simulated by keeping the surface fixed and rotating the 
ion beam correspondingly. 

In a single beam setup, a typical run of the above de- 
scribed MC simulation, starts from a flat surface, after 
t ^ 3 the first, short ripples appear, and after t ~ 6 reg- 
ular ripples of the size of the system have formed (see 
Fig.[2fa)). As fluence increases, the ripple pattern shows 
less defects (branching and deviation from expected di- 
rection). We also observe ripple motion, but will not 
report details about its characteristics in this work.^- 



III. RESULTS AND DISCUSSION 

A. Opposed ion-beams 

It has been proposed by Carter^ that IBS with two di- 
ametrically opposed ion-beams (same and a difference 
of 180° in (f)) can lead to a cancellation of instabilities 
induced by each beam. Instead, he predicted tempo- 
ral oscillations of ripple amplitudes. Furthermore, ripple 
motion should be suppressed due to the restoration of 
reflection symmetry, which would be broken by a single 
beam. 

Fig. [2fb) shows results of topographies obtained from 
MC simulations of this setup. We never observe the pre- 
dicted behavior, instead the growth of structures resem- 
bles that obtained in a single beam setup, as can be seen 
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FIG. 2: Rippled surfaces after 8ML of sputtering, (a) by a 
single ion beam, (b) by two ion beams opposed to each other. 
Arrows indicate the direction of ion beams. 
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FIG. 3: Distribution function of the asymmetry parameter e 
(see main text) for the cases of a single ion beam (upper panel) 
and two ion beams opposed to each other (lower panel). 



from comparing Figs. [2] (a) and (b). But ripples obtained 
from opposed beams appear longer and more straight 
than those from single beams and the patterns contain 
less defects. Furthermore, the shape of ripples changes 
significantly towards more symmetric slopes. In Fig. [3l 
we show a quantitative analysis of ripple slope angles ai 
and 0^2, which are defined in the inset of the upper panel. 
To measure the symmetry of the slopes, we consider the 
asymmetry parameter 



a2 — Q^i 
a2 + ai 



(1) 



Histograms of this quantity are shown in Fig. [3] for sin- 
gle beam and opposed beam setups. The average value 
of e is e ~ —0.08 for single beams (indicating that the 
steeper slope is facing the beam), and e ^ 0.01 for op- 
posed beams. We also checked, that the skewness of the 
e-distribution is reduced by using opposed beams. Thus, 
this setup may have merits in producing ripple structures 
of higher quality and order. 



B. Crossed ion beams 

Setups of multiple ion beams incident from different 
directions have been proposed by Carter in Refs IH and 
0. Vogel and Linz developed a continuum theory for a 
setup with four beams , all with the same polar angles of 
incidence and crossing with = 90°. Their theory is 
based on a damped Kuramoto-Sivashinsky equation and 
takes into account terms up to fourth order in a gradi- 
ent expansion of the erosion term. They predict square 
patterns, which can be stable in the long time regime 
under specific conditions.^ Joe et al. performed experi- 
ments with two ion-beams with = 73° and = 90° 
on Au(OOl).^^'^^ They find nanodot patterns with square 
symmetry (albeit with a rather low degree of order), if 
the fluxes of the two beams are precisely balanced. Oth- 
erwise, the structures develop into modulated ripples. 

In Fig. m upper row, we show simulation results for 
perfectly balanced ion beams (beam parameters (1) in 
Table H]). The middle row of Fig. 2] shows results for ion 
beams with different intensities incident from the same 
polar angle ((2) and (3) in Table H]) and the lower row 
shows corresponding results for ion beams of identical 
intensities, but incident from different polar angles and 
with different beam parameters ((1) and (4) in Table H]). 
White arrows in Fig. |4] indicate the projection of the ion 
beams into the x — y plane (thus they enclose A^), and 
we have chosen the geometry such that the x-axis always 
is the bisector of this angle. For balanced beams and 
< A(/) < 90°, we observe ripples with a wave vector 
parallel to the x-axis, for 90° < A^ < 180°, the ripple 
wave vector is parallel to the y-axis, but the ripple pat- 
tern shows more defects. Exactly at Acj) = 90°, square 
patterns replace the ripples. For unbalanced beams, we 
observe ripples in oblique directions, which we will dis- 
cuss below. The erosion rates parallel and perpendicular 
to the beam projection onto the x — y plane, and 
are given in Table HI The middle and lower row of Fig. |4] 
are given to emphasize that we did not observe any differ- 
ences in pattern formation due to different mechanisms 
of imbalance. All mechanisms we tried lead to the same 
patterns, if they imply the same erosion rates of linear 
theory, up to a common constant factor. 

Note, however, that in all our simulations both rates 
i^ll and u± are negative — unlike the situation in Ref. 
— indicating that the crossed beams do not mutually 
suppress their generated ripples. Thus we have nothing 
to say here about the puzzling results of Ref. UH, who find 
structures in DIBS setups, for which a straightforward 
application of continuum theory would predict smooth 
surfaces. 

Let us analyze our findings within the simple frame- 
work of linear Bradley-Harper theory. Suppose beams 



1 and 2 are characterized by erosion rates u 



1,2 



and A =||, ^. Averages of the rates over the two beams 
are denoted by ua = (^^^^ + ^^^^)/2 and we introduce 



(1) 

A 



). According to linear theory, the 
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TABLE I: Beam parameters for our DIBS setup 





e 


a 


^11 


^± 


(1) 


50° 


9.3 


-4.4 


-0.86 


(2) 


50° 


9.3 


-5.8 


-1.14 


(3) 


50° 


9.3 


-2.9 


-0.57 


(4) 


35° 


7.3 


-1.5 


-0.81 



growth of Fourier modes \h{kx^ky^t)\ ex exp(rt) is con- 
trolled by the (real) growth rate Tikx^ky)^ which is a 
quadratic form of the wave vector, i.e. V = k^z>k with a 
matrix z> of erosion rates. The growth rate of the fastest 
growing mode and its direction are obtained by deter- 
mining the largest eigenvalue of z> and the corresponding 
eigenvector. z> is easily calculated. Its matrix elements 
take on the form 

Oxx = 2(i^|| + i^x) + (^11 - ^±)cosA(/) (2) 

Oyy = 2(i^|| + i^x) - (^11 - cos A(/) 

Oxy = {Au\\ - Au±)s'mA(j), 

and Oyx = Oxy. For balanced beams, the Aua vanish 
and the erosion rate matrix becomes diagonal, indicating 
that ripples will only appear with wave vectors either 
parallel to the x-axis or parallel to the y-axis. At Acj) = 
90°, the rates in both directions become degenerate and 
square patterns will emerge, if they are stabilized by the 
nonlinear terms. As cos(A(/) + 7r/2) = — cosA(/), it is 
obvious from Eq. [2] that the regime 90° < A^ < 180° 
can be mapped to < A<p < 90° by interchanging x and 
y. This explains the main features of the upper row of 
Fig. 12 

For unbalanced beams, the orientation of ripples will 
generally depend both on A^ and the imbalances in 
growth rates Aiy\\^Au±^ but for A<p = 90° the situation 
is simpler. The cosine terms vanish, and the eigenvec- 
tors e± oc (1,±1) of the simpler matrix become inde- 
pendent of the rate imbalances. Thus, the linear theory 
predicts ripples with wavevectors inclined hy ip = 45° or 

= 135° with respect to the x-axis, i.e. parallel to one 
of the beams (from the eigenvalues it follows that the 
wave vector is parallel to the dominant beam), irrespec- 
tive of the amount or nature of the imbalance in growth 
rates, if A0 = 90°. This result is in accordance with 
the experimental finding in Ref. Il3|. The directions for 
other values of A^ are easily obtained, if the imbalance is 
known. In Fig. |4l we have indicated these directions with 
dashed white lines and find an overall satisfactory agree- 
ment of simulation results with this prediction. Fig. [5] 
shows the predicted deviations of the orientation of ripple 
wavevectors from the x-axis, i.e. the bisector of the beam 
directions for different ratios of rates, / = z^^^V^a''- 




FIG. 4: upper row: Rippled surfaces after 40ML sputtering 
by two balanced ion beams incident from = 50° and sep- 
arated in azimuthal angle by (a)A(/) = 60°, (b) A(/) = 90° 
and (c) Ac/) = 120°. middle row: same as (a)-(c) but for 
imbalanced ion beams, which differ in intensity by a factor 
of 2 (/ = 1/2) corresponding to beams (2) and (3) in Ta- 
ble [H lower row: same as middle row, here the imbalance is 
generated by different beam parameters, (1) and (4) in Ta- 
ble ID Arrows indicate the directions of ion-beams. In middle 
and lower row bigger arrows correspond to the dominant ion 
beam. Dashed lines indicate the directions predicted for the 
wavevectors of ripples by linear theory (see main text). 



C. Sputtering of rippled surfaces 

Since the possibilities of simultaneous sputtering by 
multiple beams are often limited in experiments, sequen- 
tial sputtering by a single beam from different directions 
seems to be a more promising setup. Vogel and Linz 
proposed SIBS as a general substitute for multi-beam 
sputtering.^ Note that in SIBS setups, a precise balance 
between fluxes of multiple beams, which may be diffi- 
cult to achieve in DIBS, can be adjusted by tuning the 
exposure time in each direction. 

Joe et al.^^ have recently performed experiments on 
Au(OOl) using an ion beam (incident from = 72°), for 
which ripples with wave vectors perpendicular to the pro- 
jection of the beam direction into the initial surface plane 
are formed. After rotating the target by 90° (keeping 6 
fixed), the initially formed ripples are very rapidly de- 
stroyed and new ripples build up in the correspondingly 
rotated direction, but the authors could not get patterns 
corresponding to a superposition of two generations of 
ripples as expected from the predictions of continuum 
theory. 

We have simulated SIBS with ion beams incident from 
= 50° and a rotation step of Acj) = 90° after 9 monolay- 
ers of erosion. Initially, ripples appear with wave vectors 
parallel to the direction of the projected ion beam into 
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FIG. 5: Orientation of fastest growing mode for different val- 
ues of A(/). ip denotes the angle of the ripple wave vector with 
the X-axis. Different lines correspond to different imbalances 
f (see main text). 



the x-y plane (see Fig. [6fa)). After the rotation step, 
a correspondingly rotated ripple pattern builds up as 
shown in Fig. [6fc) . In a narrow time window, shortly af- 
ter the rotation step (shown in Fig.[6l^b), at t = 10.8ML) 
a superposition of ripples of both orientations is observed. 
Fig. injd) shows the structure factor of the height profile, 
which depicts the degree of order. 

Beyond o:^ IML of erosion after the rotation step the 
time evolution of surface roughness (shown in Fig. [7j) 
differs strongly from what one expects in linear contin- 
uum theory. There, the height of prestructured ripples, 
having already reached saturation, would neither grow 
further, nor will it decrease. Ripples generated by the 
new beam direction will grow exponentially and thereby 
catch up the height of the prestructure. The different 
scenario we observe in the simulations is shown in the 
inset of Fig. [71 After IML, during which the system 
follows the scenario predicted by Bradley-Harper theory, 
the prestructure collapses very rapidly and is flattened, 
leading to a transiently decreasing roughness. We define 
a structural relaxation time Tq as the length of the inter- 
val from the rotation step through the transient decrease 
of roughness to the point, where roughness has retained 
its value immediately before the rotation step. Roughly 
this time interval contains all processes necessary to ro- 
tate the ripple pattern from the previous to the current 
rotation step. In our case Tq 2.5. Although Tq can- 
not be obtained from linear rate theory our observations 
need not be in contradiction with continuum theories, as 
the prestructure constitutes an initial condition beyond 
the range of validity of linear theory. Let us reemphasize 
that the roughness decreases transiently, despite the fact 
that both linear erosion rates and iy± (see Table (|T|)) 
indicate an unstable growth of fluctuations in x and y 
direction. A superposed square pattern is observable, if 
the height of the growing new ripples and the shrinking 




FIG. 6: Snapshots of surface profile during SIBS. At t = 9, 
a A(/) = 90° rotation step occurs. Shots are at (a) t = 9, (b) 
t = 10.8 (inset: zoomed in) and (c) t = 18. Arrows indicate 
the direction of ion beams. In panel (d) the structure factor 
\h{kxj ky)\'^ of the profile of panel (b) is shown. 




FIG. 7: Roughness of surface patterned by SIBS vs. time. 
The surface is rotated by 90° at t = 9. The inset shows the 
temporal evolution of the maximum of the structure factor 
Smax = max(|/i(S)|^) for wavevectors in || and in ± direction. 



old ones become comparable. For our parameters this 
takes place at ~ 1.5 ML after rotation in a narrow time 
window. In the experiments by Joe et alM values of 
growth rates are = 14.2 and Uy = —0.9. Since Ux > 
and Uy < 0^ linear rate theory predicts that the old gen- 
eration of ripples will decay immediately after a rotation 
step, which may cause an even faster collapse of the pre- 
structure and make the time window of existence of a 
superposed structure unobservably small. 
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D. Sputtering of continuously rotating sample 

Sample rotation during IBS is applied for various rea- 
sons. One motivation is to achieve suppression of pattern 
formation in SIMS and AES depth profiling/'^ There are 
numerous reports that RIBS can suppress surface rough- 
ening and enhance the resolution of depth profiling (see 
Ref. i21j and its references). Although this method is fre- 
quently used, there are very few systematic studies of 
the effects of different parameters, in particular of the 
angular velocity of rotation. Previous discrete simula- 
tions could not successfully explain the observed strong 
suppression of roughening. Recently, IBS has become a 
popular method for smooth etching of metallic surfaces. — 
Here too, sample rotation has been proposed as a prac- 
tical measure^^ to prevent nano-scale roughening. A dif- 
ferent motivation to use RIBS is the control of pattern 
formation. Frost et al. have found that off-normal IBS 
with sample rotation may lead to formation of hexago- 
nal, close packed quantum dots i^^i^^ This is attributed to 
a restoration of axial rotation symmetry with respect to 
the average surface normal, which is broken by off-normal 
incidence of a single ion beam at fixed azimuthal angle. 
Dot formation in rotated, off-normal IBS has been found 
in the framework of continuum theories^^^^^^^, assuming 
a fiux of incoming ions, which is distributed evenly over 
all azimuthal angles. This assumption corresponds to 
the limit of high rotation frequencies. Dots also appear 
in MC simulations performed in the high rotation fre- 
quency Hmit for a wide range of parameters. 

In the present work, we focus on the systematic depen- 
dence of height fiuctuations on the rotation frequency, 
irrespective of the random or deterministic nature of 
these fiuctuations. Therefore, we will study the rough- 
ness w = {{h — (h))'^) ^ averaged over the sample and an 
ensemble, as function of time and rotation frequency uj. 
Reported rotation frequencies cover a range from 0.1 to 
15 rpm for different fiuxes and different types of ions and 
materials jii^»2£i2^i2£ There is a predictions of the scaling 
of height with cj, which, in Ref. 131|, is given in the form 



h{k^ t) (X exp( 



- ^11 ) 



sin(2a;t)). 



(3) 



It is based on the original linear Bradley-Harper theory 
in rotating coordinate systems. Furthermore, Cui et al. 
reported that ripples do not form for angular frequencies 
greater than 0.1 rpm.^^ (for a fiux of 3.5 x 10^^ ions cm~^ 
s-^ of 300 e V Ar+ on GaN substrate). 

In section IIII C\ we have defined a characteristic re- 
sponse time To(A(/)) of the IBS generated structures to 
sudden changes of the beam direction. We propose 
that this time scale is also of relevance for RIBS, as sur- 
face structures might follow rotation frequencies much 
smaller than ujq = A^/Tq adiabatically. On the other 
hand u Uq might correspond to the high frequency 
limit. We performed simulations with different rota- 
tion frequencies varying in the range 0.05 • • • 50co'o, with 
ujo 36° per eroded monolayer taken from our SIBS sim- 
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FIG. 8: Time evolution of surface roughness during RIBS for 
different rotation frequencies. 



ulations. For low frequencies, ripples form and rotate in 
synchrony with the beam direction. For high frequen- 
cies, (cj = oo is included as it corresponds to random 
azimuthal directions of incoming ions chosen from a flat 
distribution), cellular structures of growing size are ob- 
served. These flndings are in agreement with predictions 
in Ref. 9. A mixture of short ripples and cellular struc- 
tures appears at intermediate uj. Somewhat surprisingly, 
the roughness of the surface is not a monotonous func- 
tion of rotation frequency. In Fig. [8] roughness is shown 
as a function of time for different uj. For low (includ- 
ing UJ = 0) and high frequencies, it grows monotonically, 
approaching approximately equal growth rates beyond 
t ^ 30, independent of uj. For intermediate cj, the rough- 
ness shows oscillations (with frequencies Q ~ 4co'), and 
its growth rate is strongly suppressed. For uj = 3.5 ujq 
growth rate reaches a minimum. In Fig. [9] the roughness, 
averaged over a period T = 2i\ juj^ is shown against cj, for 
times up to t = 37 . For increasing the minimum at 
UJ 3.5 cjo becomes more and more pronounced. Thus 
our simulations predict an optimal rotation frequency, 
if preparations aim at smooth surfaces. This frequency 
also marks the frequency scale beyond which ripples do 
no longer appear. The very rapid crossover to non-ripple 
structures is in accordance with flndings of Ref. [30|- 

To compare our results with the prediction of Eq. [3] 
about the scaling of height with cj, we studied S = 
|/i(/c,t)p for a flxed value of wave vector k. It grows 
rapidly and then oscillates with frequency 2uj around a 
saturation value with an cj-dependent amplitude c, as 
shown in the inset of Fig. [TOl The oscillatory behavior 
with frequency 2uj is also present in Eq. [3l The main part 
of the flgure depicts the decrease of the oscillation ampli- 
tude with increasing a; in a double logarithmic plot. The 
fltted line has a slope of — 1.05 ± 0.05, which is very close 
to the 1/uj behavior suggested by Eq. [3l Note, however, 
that this equation was derived from linear theory and 
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FIG. 9: Roughness against angular frequency uj at different 
times. 




FIG. 10: The logarithm of the oscillation amplitude of struc- 
ture factor for a given /c, vs. loguj. Inset: Time evolution of 
structure factor for a given /c, for uj — 2uj{) displays the oscil- 
lations. The amplitude of these oscillations is shown in the 
main figure. 



has to be multiplied by an exponential growth factor, 
whereas our result applies to the non-linear saturation 
regime. Due to the rapid initial growth, the regime of 
validity of Eq. [3] remained unobservable. 

We also measured the total amount of eroded material, 
AM up to t = 40 for different rotation frequencies. This 
integrated yield decreases with increasing uj and displays 
a clear distinction between a low- and a high-frequency 
regime in a semi-log plot, shown in Fig.[TTJ The interpo- 
lated crossover frequency between these regimes is very 
close to S.ScJo, the frequency, which minimizes the to- 
tal roughness (see Fig. [9]). This is in accordance with 
our findings (see Fig. [7]) that more pronounced rippled 
prestructures (developing at low rotation frequencies at 
an angle uot) decay more rapidly (at t + dt) than the 




FIG. 11: Total amount of eroded material (integrated yield) 
after 40ML of sputtering of rotating samples vs. rotation 
frequency. The arrow is at :^ 3.5cJo, the frequency which 
minimize the total roughness. 



smaller and less regular height fluctuations, which ap- 
pear at higher rotation frequencies. 



IV. CONCLUSION AND OUTLOOK 

We studied dual ion beam sputtering (DIBS), sequen- 
tial ion beam sputtering (SIBS) and rotating ion beam 
sputtering (RIBS) by a kinetic MC simulation technique, 
which combines erosion events due to single ions and 
surface diffusion. For a DIBS setup with two diamet- 
rically opposed beams, we did not confirm predictions by 
Carter^, but rather found non- moving ripples with orien- 
tations as in a single ion beam setup. The ripples have 
a higher degree of order and more symmetrical slopes as 
compared to those created by single beam sputtering. For 
DIBS setups with crossed ion-beams, we find ripple pat- 
terns for crossing angles A(j) ^ 90° and square patterns 
for crossing of balanced beams at exactly right angle. 
The ripple orientations follow the predictions from lin- 
ear Bradley Harper theory. Any kind of beam-imbalance 
leads to ripple patterns oriented according to the domi- 
nant beam. This is in accordance with the experimental 
observation in Ref. 13. For SIBS setups, we found a 
vary rapid destruction of the ripple prestructure of the 
previous rotation step, which cannot be explained by lin- 
ear Bradley-Harper theory. The flattening of the pre- 
structure leads to a transient decrease in total roughness. 
Only within a very short time window, the growing new 
generation of ripples and the shrinking old ones lead to 
a superposed square pattern. Thus we could not con- 
firm propositions to use SIBS as a universal substitute 
for complicated multi-beam setups. The rapid destruc- 
tion of the prestructure is in accordance with findings of 
Ref. 13. For RIBS setups we observed a non-monotonic 
dependence of roughness upon rotation frequency. At 
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a frequency scale set by the structural relaxation time 
of prestructures, which can be observed in SIBS simula- 
tions, an increasingly pronounced minimum of roughness 
occurs with increasing time. We found that the struc- 
ture factor at fixed wavevector rapidly approaches sta- 
tionary oscillations around a saturation value with os- 
cillation amplitudes inversely proportional to frequency. 
This behavior was also predicted from linear theory, but 



seems to have a much broader range of validity. 
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